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The occurrence of energetic and dynamical instabilities in a Bose-Einstein condensate moving in 
a one- dimensional (ID) optical lattice is analyzed by means of the Gross-Pitaevskii theory. Re- 
sults of full 3D calculations are compared with those of an effective ID model, the nonpolynomial 
Schrodinger equation, pointing out the role played by transverse degrees of freedom. The insta- 
bility thresholds are shown to be scarcely affected by transverse excitations, so that they can be 
accurately predicted by effective ID models. Conversely, transverse excitations turn out to be im- 
portant in characterizing the stability diagram and the occurrence of a complex radial dynamics 
above the threshold for dynamical instability. This analysis provides a realistic framework to discuss 
the dissipative dynamics observed in recent experiments. 
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I. INTRODUCTION 

The occurrence of energetic and dynamical instabilities 
in a Bose-Einstein condensate (BEC) that moves in a pe- 
riodic (optical) potential is an interesting problem from 
the conceptual viewpoint, since it involves basic prop- 
erties of superfluids. Experimentally the dissipative dy- 
namics of BECs moving through optical lattices has been 
investigated both in the weak 0, and tight-binding 
regimes 0, , and a number of papers have recently been 
published in connection with these and similar experi- 
ments HHHEHESmHIil. 

From the theoretical side, systematic investigations 
of the stability regimes have been so far presented for 
strictly one-dimensional (ID) systems. For deep optical 
lattices (tight-binding limit) it has been argued that the 
major mechanism responsible for the superflow break- 
down of a moving BEC is the onset of a dynamical 
(modulational) instability 00- This behavior has also 
been confirmed by direct integration of the 3D Gross- 
Pitaevskii (GP) equation in the same regime [lll[T21 |. and 
well reproduces the experimental observations ]4$. Con- 
versely, the opposite limit of shallow lattices has been the 
object of a stimulating debate and it has been suggested 
that the origin of the observed instability could have an 
energetic or dynamic character 0,0,0]. 
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Here we present a general discussion of the energetic 
and dynamical instabilities within the GP theory. We 
make use of a 3D description that includes the radial de- 
grees of freedom in order to provide a reliable framework 
for a comparison with current experiments. In Sec. II, 
we briefly summarize the formalism of the standard linear 
analysis within the GP theory. Then, in Sees. Ill and IV, 
we discuss in detail the case of a cylindrical condensate in 
a ID optical lattice, for which one can use Bloch functions 
and rigorously define the concept of quasimomentum as- 
sociated with the motion of the condensate in the lattice. 
For any given quasimomentum of the condensate, we cal- 
culate the excitation spectrum and the stability diagram 
by solving the equations of the linearized GP theory . 
The same quantities are also calculated by means of an 
effective ID model, known as nonpolynomial Schrodinger 
equation (NPSE) |15|], which includes the transverse di- 
rection through a Gaussian ansatz for the radial shape of 
the order parameter, with a z- and i-dependent width. 
With respect to strictly ID models, which rely on a suit- 
able renormalization of the mean- field coupling constant, 
the NPSE has the advantage that the true 3D coupling 
constant g = 47th 2 a/m can be used, thus allowing for a 
direct comparison with experiments. It turns out that 
this model gives accurate predictions for the instability 
thresholds, the latter being mainly determined by the 
dispersion of the lowest branch of excitations, with no 
radial nodes. Higher branches, with one or more radial 
nodes, that are included in the GP theory but not in the 
NPSE, are shown to be important in characterizing the 
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number and the type of excitations that become unstable 
above the instability threshold. They also produce a sig- 
nificant change in the shape of the dynamically unstable 
regions above threshold. 

Finally, in Sec. V, we perform a GP simulation for the 
elongated condensate of the experiments of Ref. 0] . The 
results of the simulation are analyzed by using the spec- 
tra and the phase diagrams of the cylindrical condensate 
of Sees. Ill and IV. This analysis provides convincing 
evidence that the breakdown of superfluid flow observed 
in Ref. Q is associated with the onset of a dynamical 
instability, triggered by the resonant coupling of Bogoli- 
ubov phonon and antiphonon modes. The subsequent 
dynamics is shown to strongly involve radial motions. 



II. GP EQUATION AND LINEAR ANALYSIS 

The GP equation for the order parameter ip oi N con- 
densed atoms of mass m and scattering length a is [l4| 
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where g = Airh a/m. We will consider the case of con- 
densates confined in a harmonic potential superimposed 
to a ID periodic optical lattice, V(x) = T4o(x) + Vl(z). 

Let ipo(x,t) — </>o(x) exp(— ifit/h) be a solution of the 
stationary GP equation, having chemical potential /i, 
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and normalized according to Jd 3 x |0o( x )| 2 = N. If the 
state 4>q is a local minimum of the energy functional 
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the system is energetically stable. Energetic (Landau) 
instability sets in when this condition is no longer true, 
and there are some directions in the ip space along which 
it is possible to lower the energy. This happens, for in- 
stance, when a uniform condensate flows with velocity v 
in presence of a static external potential; if v is larger 
than the velocity of Bogoliubov's sound, then the sys- 
tem can lower its energy by emitting phonons (Landau 
criterion of superfluidity 16]). Similar conditions can be 
found for an inhomogeneous condensate moving in an 
optical lattice above a certain critical velocity Ve- 

The condition for energetic stability can be formally 
derived by considering small deviations from 0o in the 
form 
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and expanding the energy functional up to the quadratic 
terms 
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The onset of energetic (Landau) instability is thus sig- 
naled by the appearance of negative eigenvalues in the 
spectrum of the operator M. 

A second type of instability is the dynamical (modu- 
lational) instability. This occurs when the frequency of 
some modes in the excitation spectrum has a nonzero 
imaginary part. Then the occupation of these modes 
grows exponentially in time and rapidly drives the sys- 
tem away from the steady state. The conditions for the 
onset of this instability can be obtained again by means 
of the expansion (0} . When inserted into the GP equa- 
tion (JTJ it yields the Bogoliubov equations 
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whose character is now determined by the operator a z M, 
with 
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Energetic and dynamical instabilities are therefore re- 
lated to the spectral properties of the operator M and 
a z M 0,0. This has two important consequences: (i) 
when the system is energetically stable it is also dynami- 
cally stable; (ii) the energetic instability, which is related 
to the spectrum of M, cannot be revealed by a direct inte- 
gration of the time-dependent GP equation (JTJ , which is 
dissipationless and is governed by the spectrum of a z M. 



III. SPECTRUM OF A CYLINDRICAL 
CONDENSATE IN A LATTICE 

In order to perform explicit calculations, let us con- 
sider the case of an infinite cylindrical condensate, which 
is radially confined by the harmonic potential Vh (r) = 
(l/2)ma-' 2 r 2 , and is subject to the periodic potential 
Vl(z) = sE R cos 2 (qsz). The Bragg wave vector qs = 
ir/d is determined by the lattice spacing d. The quantity 
En = h 2 qj 3 /2m is the recoil energy of an atom absorbing 
one lattice photon and s is a dimensionless parameter fix- 
ing the lattice intensity. Such a cylindrical condensate is 
an accurate representation of an elongated axially sym- 
metric condensate, when the weak axial harmonic con- 
finement can be neglected. The major advantage of this 
geometry is that one can exploit the periodicity of the 
system by expanding the order parameter on a basis of 
Bloch functions, with the same period of the lattice. 



A. Linearized GP equation 

Within the GP theory, the starting point of the lin- 
ear stability analysis for the cylindrical condensate is the 
expression 

^(r, z, t) = e-WptjP* [<t> p0 (r, z) + %(r, z, t)} . (10) 

The Bloch wave vector p (also called quasimomentum 
of the condensate) is associated with the velocity of the 
condensate in the lattice and is restricted to the first 
Brillouin zone 0, @ • For a given p, the function <p p0 is 
the solution of the stationary GP equation with energy 
E p , while the excited part of the order parameter can be 
written as 

8<t> p {r, Z ,t) = Y.( u p^ z y {qZ ~ UJpq ^ ) 

1,3 

+v; qj (r 7 z)e-^ z ~ UJ P^) , (11) 

where q is the Bloch wave vector (or quasimomentum) of 
the excitations, and j represents all possible values of the 
two quantum numbers a = 1, 2, 3, . . . and v = 0, 1, 2, . . .. 
The former is Bloch band index and the latter is the 
number of radial nodes in the Bogoliubov quasiparticle 
amplitudes u and v. Along z all the functions <p P Q, u pq j 
and v pq j are periodic with the same period of the lattice. 

Inserting the above expressions into the GP equation 
(JTJ, one gets the Bogoliubov equations 
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The quasiparticle amplitudes u and v satisfy the relation 
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Note that, for each quasiparticle obeying this orthonor- 
mality conditions, the Bogoliubov equations admit an 
antiquasiparticle with norm — 1, which corres pon ds to 
replacing the amplitudes (u, v) with (v*,u*) |7lll7j. 

The above equations can be solved numerically. We 
choose the relevant parameters such to simulate typical 
elongated condensates of the experiments of Ref. . The 
linear atomic density in the z direction, averaged over the 
lattice length d, is taken to be equal to the linear density 
of the actual 3D condensate, close to the trap center and 
in the absence of the lattice. Here the chemical poten- 
tial for p = is taken to be fio = 7.2 huj r . We also use 
2-k qg l — 800 nm and u> r = 2tt x 90 Hz. For each p, the 




FIG. 1: Excitation spectrum of a cylindrical condensate at 
rest in the optical lattice (quasimomentum of the condensate 
p — 0) as a function of the quasimomentum q of the excita- 
tions, for lattice intensity s = 5. The first two Bloch bands are 
shown and for each band we plot the first 20 radial branches. 
The lowest branch in each band (thickest line) corresponds to 
axial excitations with no radial nodes. In the inset we show 
a magnification of the low energy spectrum. 



order parameter (f> p o(r,z) and its energy E p are calcu- 
lated by solving the stationary GP equation by means of 
iterative methods. Then the Bogoliubov equations 1|12|) 
and l|13|) are transformed into matrix equations by pro- 
jecting the functions u and v on the Bessel-Fourier basis 
\l, n) = J (a Qn r/r max ) exp(i2q B lz), where a Qn are zeros 
of the Bessel functions Jn(r) and r max is the radius of 
the computational box [l^, [Hj • The resulting equations 
are not diagonal over the Fourier index I, and involve the 
numerical diagonalization of a 2 x N r x N z rank matrix. 

The p = case corresponds to the condensate at rest. 
In Fig. n w e show a typical spectrum, for s — 5. The first 
two Bloch bands (a — I and 2) are shown. For each band, 
we plot the first 20 radial branches (v = 0,1,2, ... , 19). 
The thickest lines are the dispersion laws of the v — ex- 
citations. A magnification of the low-energy branches is 
shown in the inset. The lowest branch corresponds to the 
dispersion law of v = axial phonons. It starts linearly 
at q — > and its slope is the Bogoliubov sound velocity c. 
We calculate the slope for different values of the lattice 
intensity s and we plot the results in Fig. [21 (points). The 
dashed line is the analytic prediction c = (w*k) -1 / 2 of 
Ref. [2(j , where the effective mass m* and the compress- 
ibility k are defined as m* — \im p ^ (d 2 E p /dp 2 )~ 1 and 
= n(dp/dn) (see also the discussion in Ref. 13J). 
Here we calculate the effective mass from the numerical 
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FIG. 2: Bogoliubov sound velocity c as a function of the lat- 
tice intensity s. Points: slope of the lowest phononic branch 
in the Bogoliubov spectrum of the cylindrical condensate. 
Dashed line: analytic prediction c = (Km*) -1 ' 2 (see text). 



values of the energy E p . The compressibility is taken to 
be the one of a uniform gas in the same ID optical lat- 
tice and with density equal to the average density of the 
cylindrical condensate (for the radial average we use the 
Thomas- Fermi value n — n(Q)/2 where n(0) is the cen- 
tral density). The agreement is remarkable. The value 
at s = is the sound velocity c = \qn(0)/ (2m)] 1 / 2 in a 
cylindrical condensate without lattice |2l|. Finally, we 
note that the v = 1 branch in the inset of Fig. 2] exactly 
starts at u> = 2uj r for q = 0, where it corresponds to the 
radial breathing mode of the cylindrical condensate. 

At p 0, when the condensate moves in the lattice, the 
whole spectrum changes. In Fig. [3] we show a sequence of 
six spectra for s = 5 and p/qb = 0, 0.25, 0.5, 0.55, 0.75, 1. 
The real part of the frequencies io pq .j is plotted. At p = 
quasiparticles and antiquasiparticles have the same en- 
ergy, since the system is symmetric under z — > —z in- 
version. For a given q, they correspond to excitations 
propagating in opposite directions (opposite sign of lu). 
In Fig. |21 the spectrum of antiquasiparticles, at p = 0, is 
plotted in the negative ui plane. 

By increasing p, the slope of the v — phononic branch 
increases for phonons which propagate forwards on top of 
the traveling condensate. Vice versa it decreases for those 
propagating backwards and their dispersion law eventu- 
ally crosses the u> = axis and changes sign. By fur- 
ther increasing p, one notices that the frequency of both 
phonons and antiphonons approaches zero at the bound- 
aries of the Brillouin zone. At a critical quasimomentum 
Pd {pd — 0.525 qs in our case), the real part of the 
v = mode frequency vanishes at q = iqs- Phonon 
and antiphonon pairs thus exhibit a resonance coupling 
by first-order Bragg scattering, giving rise to a dynami- 
cal instability. For ID condensates this process has been 
already discussed in Ref. Q. The coupled pairs of ex- 
citations, with \q\ = qB = n/d, start having imaginary 
frequency. The system can therefore develop a macro- 
scopic standing wave, whose wavelength is twice the lat- 
tice spacing, which implies an out-of-phase oscillation of 




0.5 0.5 

q/q B 



FIG. 3: Real part of the excitation spectrum of a cylin- 
drical condensate in a lattice with s = 5 and for differ- 
ent values of the condensate quasimomentum p (p/qB = 
0,0.25,0.5,0.55,0.75, 1), calculated from Eqs. iPt - ljT^ . 



the condensate population in adjacent lattice sites. This 
is also consistent with the results of the ID simulations 
of Ref. @. 



For p > po the coupling between v — phonons and 
antiphonons extends from zone boundary to lower values 
of q. Moreover, modes with radial nodes (y > 0) also 
couple, first at q = qB and then down to lower q's. A 
conjugate pair of complex frequencies appears each time 
a resonant coupling occurs between a pair of quasiparti- 
cle and antiquasiparticle that are degenerate prior to the 
coupling. An example is given in Fig. 01 where we plot the 
real and imaginary parts of the eigenfrequencies obtained 
by solving Eas. i|12fl and (|13f) in the case p/qB — 0.55. Fi- 
nally, as one can see in the last two frames of Fig. |3| a 
complicate sequence of no-crossing patterns also develop 
between modes belonging to the quasiparticle spectrum 
or between those belonging to the antiquasiparticle spec- 
trum. These no-crossings do not produce any dynamical 
instability, but they contribute to mix up several radial 
branches via hybridization processes, which are not in- 
cluded in ID models. 
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FIG. 4: Real (top) and imaginary (bottom) parts of the 
frequencies obtained from Eqs. 112H - I13|I for s — 5 and 
p/qB = 0.55. The vertical (dotted) lines correspond to the 
points where a conjugate pair of complex frequencies appears. 



B. Nonpolynomial Schrodinger equation 



The description of the system can be simplified 
by using an effective ID model, the nonpolynomial 
Schrodinger equation (NPSE) which partially in- 

cludes also the radial to axial coupling, and has been 
shown to pro vide a realistic description in several situa- 
tions [Hd!|22. This model is obtained from the GP 
equation by means of a factorization of the order param- 
eter in the product of a Gaussian radial component of 
z- and ^-dependent width, a(z,t), and of an axial wave 
function tp(z,t) that satisfies the differential equation 
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coupled with an algebraic equation for the radial width 
a[z, t) = a r ^/l + 2aN\^(z,t)\ 2 . (18) 



Here a r — y/KJ mio r is the radial oscillator length. By 
expanding ?p(z,t) around the stationary wave function 
4>o( z ) and using Eqs. (|17fl and l|18l) . it is straightforward 
to obtain the Bogoliubov-like equations 
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The solution of these equations is obtained by a numeri- 
cal matrix diagonalization. Compared to the solution of 
the Bogoliubov equations i|12|) and (|13fl . the calculation 
of the spectrum from Eqs. I|19f) is much faster. 

Owing to the partial coupling between axial and radial 
degrees of freedom, the NPSE provides a more accurate 
description of the actual 3D condensates with respect to 
strictly ID models. An important advantage consists in 
the fact that, differently from strictly ID models, the 
NPSE does not require any renormalization procedure 
for the coupling constant g, thus offering the possibility 
to make a direct comparison with experiments, as well as 
with GP calculations for 3D condensates. In this respect 
one has to remark also that, since the radial degrees of 
freedom are included through a single function cr(z,t), 
the NPSE only accounts for v — radial excitations. The 
comparison between the predictions of the NPSE and 
those of the GP equation for the cylindrical condensate 
is thus particularly helpful to point out the role of the 
v > radial modes. 

In Fig. [S]we plot the real (solid lines) and imaginary 
part (dashed lines) of the excitation frequencies obtained 
with the NPSE for the same s and p of Fig. [31 The two 
branches correspond to the v = phonon and antiphonon 
dispersions and their behavior is shown to closely follow 
that of the lowest v = branches of the GP case. In 
particular, the critical value p = pd at which the two 
modes couple at q = qB, giving rise to dynamical insta- 
bility, coincides with the GP result (pu = 0.525 qs)- The 
qualitative behavior is very similar also to that found in 
the ID analysis of Ref. 0. 



IV. STABILITY DIAGRAMS 

By repeating the calculations of the Bogoliubov fre- 
quencies (eigenvalues of the matrix a z M) for different 
values of the condensate quasimomentum p, one can draw 
the dynamically unstable regions in the p-q plane. In the 
same way, one can find the spectrum of the matrix M and 
draw the regions of energetic instability. The calculations 
can be performed both with GP theory and NPSE. Typ- 
ical results are shown in Fig. for s — 1,5, and 10. The 
white area corresponds to a stable condensate; the light 
shaded area to a condensate which is energetically unsta- 
ble (in the presence of dissipative processes, the energy 
can be lowered by emitting phonons of quasimomentum 
q which lie in the shaded range); finally the dark shaded 
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FIG. 5: Excitation spectrum obtained from the NPSE for a 
cylindrical condensate in a lattice with s — 5 and for differ- 
ent values of the quasimomentum p, as in Fig. |3] (p/qa = 
0,0.25,0.5,0.55,0.75,1). 



area is the region where the system is both energetically 
and dynamically unstable. The thresholds for energetic 
and dynamical instability correspond to the lowest values 
of p at which the instability occurs. Let us call them pe 
and pd, respectively. 

Two main features emerge from the comparison of GP 
with NPSE: (i) the left border of the energetically and 
dynamically unstable regions are very similar in the two 
cases, which also implies a a good agreement for the criti- 
cal quasimomcnta pe and p^', (h) for small values of s the 
shape of the dynamically unstable region above threshold 
looks rather different. 

The first point is better visualized in Fig. \7\ where we 
show the results for pe and po as a function of s. The 
predictions of the two models arc almost indistinguish- 
able for both the energetic and dynamical instabilities 
in the whole range of s here considered. These results 
confirm that the NPSE can be used as an efficient model 
to quantitatively identify the onset of instability for a 
wide range of s, relevant for available experiments |l2j . 
The reason for its accuracy in predicting the thresholds 
is simply that the critical quasimomenta pe and po are 
determined by the behavior of the v — branch of ex- 
citations with no radial nodes. A comparison between 
Figs. 01 and shows that these modes are indeed very 
accurately described by the NPSE up to the instability 
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FIG. 6: Stability diagrams obtained from the GP spectra 
(left) and NPSE spectra (right) as a function of the quasimo- 
menta p and q of the condensate and of the quasiparticles, 
respectively, and for s = 1,5, 10. Shaded areas represent the 
regions where the system is dynamically (dark) and/or ener- 
getically (both light and dark) unstable. Dashed lines corre- 
spond to the modes having the largest imaginary part (most 
unstable modes). 



threshold. 

A closer inspection of the stability diagrams shows a 
difference in the shape and size of the dynamically unsta- 
ble region at small s. The one predicted by NPSE does 
not fill the whole region above threshold. This agrees 
with previous ID models 0,111], but disagrees with the 
predictions of 3D GP equation. This disagreement can 
be easily understood by looking at the role played by 
the transverse degrees of freedom. In the NPSE, which 
only accounts for v — modes, nonvanishing complex 
frequencies appear when phonons and antiphonons cou- 
ple, collapsing onto a single dispersion law. For a given 
p, this provides a range of q for unstable modes, which 
starts first at qs and then lowers towards q — > (see 
Fig. El). However, at high enough p, the dispersion of the 
phonon-antiphonon pair splits again into two separate 
curves and this gives a window of stable modes from a 
certain q (at fixed p) up to zone boundary. This window 
is the dynamically stable region on the top-right corner of 
the NPSE stability diagrams in Fig. El Now, by looking 
at the GP spectra in Fig. |3 one can easily understand 
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FIG. 7: Thresholds for energetic (circles) and dynamical 
(squares) instability as a function of the lattice intensity s. 
The results of GP (filled symbols) and NPSE (empty sym- 
bols) are almost indistinguishable. Shaded areas represent 
the regions where the system is dynamically (dark) and/or 
energetically (both light and dark) unstable, according to the 
GP. 

why these stable regions are absent in the GP stability 
diagrams. The v — modes in this case are mixed up 
with several v > branches and these radial modes can 
be unstable also in the range of q where the v = modes 
are stable. It is worth stressing that the occurrence of 
complex frequencies for such radial excitations also im- 
plies that, once the system is brought into the unstable 
region above po, it can easily develop macroscopic mo- 
tions with a nontrivial radial dynamics. 

Finally, it is interesting to compare the predictions of 
GP and NPSE for the growth rate, |Im(w)|, of the most 
dynamically unstable modes. The results are shown in 
Fig.|SJ At small s, the agreement between the two calcu- 
lations (solid and empty symbols) is rather good. By in- 
creasing s the agreement for the absolute values of growth 
rates becomes less and less quantitative. However, for 
each value of s, the shape of the GP and NPSE curves 
is still quite similar. This result is interesting in view 
of the comparison with experiments, where the shapes of 
these curves can be more relevant than their absolute val- 
ues. In the experiments, in fact, it is much more likely to 
measure loss rates (or lifetimes) on large timescales, and 
estimate their relative increase as a function of p, rather 
than the absolute values of the growth rates in the small 
amplitude linear regime (short times). 

V. AN EXPERIMENT REVISITED 

In this section we compare the predictions of the lin- 
ear analysis discussed so far with the direct solutions of 
the time-dependent Gross-Pitaevskii equation . As an 
example we consider the case of the LENS experiment 
presented in Ref. which has been the object of a 
stimulating debate on the origin of the observed insta- 
bility Indeed two interesting questions have been 
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FIG. 8: Growth rates of the most dynamically unstable modes 
for different values of the lattice intensity s (s = 0.1 triangles, 
s = 1 squares, s = 5 circles) obtained with both GP and 
NPSE (filled and empty symbols, respectively). 



raised: (i) whether the GP equation is suitable to de- 
scribe the phenomenology observed in the experiment, 
and (ii) whether the observed breakdown of superfluid 
flow comes from energetic and/or dynamical instability. 

To model the LENS experiment pj, we consider an 
elongated 87 Rb condensate with N = 3 x 10 5 atoms con- 
fined in an axially symmetric harmonic trap of frequen- 
cies u r — 2-7T x 90 Hz, uj z = 2tt x 8.7 Hz. The optical 
lattice along z is produced with laser beams of wave- 
length A = 2ir/qB — 795 nm and intensity s = 1.59. The 
system is prepared in the ground state of the combined 
harmonic+periodic potential and then is let evolve after 
a sudden displacement Az of the harmonic trapping. 

If one considers not too large displacements the evo- 
lution of the condensate wave function can be safely de- 
scribed by a state of well defined quasimomentum (before 
instabilities occur) [l^. Using the formalism of Sees. Ill 
and IV we first calculate the stability diagram predicted 
by the linearized GP theory for an infinite cylinder having 
the same linear density as the one of the actual 3D elon- 
gated condensate near its center. The diagram is plotted 
in the upper part of Fig. [5] In the lower part we show 
the center-of-mass velocity, defined as v cm = dE p /dp, ob- 
tained by solving the stationary GP equation (solid line). 
The two horizontal dotted lines indicate the critical veloc- 
ities for the onset of energetic (ve) and dynamical (i>d) 
instability. 

Then we numerically integrate the full 3D GP equa- 
tion, using the experimental parameters. The evolution 
of the condensate and of its center-of-mass velocity, af- 
ter an initial lateral displacement of 60 /xm, is shown in 
Fig. As the figure shows, the condensate first acceler- 
ates towards the new trap center (to the left in the insets) 
by keeping its shape. The overall phase of the order pa- 
rameter is also preserved, except for the addition of the 
extra phase associated with the translational motion. At 
a critical velocity the coherence is suddenly lost and the 
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FIG. 9: Top: GP stability diagram for a cylindrical conden- 
sate which approximates the experimental elongated conden- 
sate of Ref. [J in its central region. Light and dark shaded 
areas are unstable regions, as in Fig. Bottom: center-of- 
mass velocity as a function of the condensate quasimomentum 
(vb = TiqB/m). The critical velocities for the energetic (ve) 
and dynamical (vd) instabilities of the cylindrical condensate 
are shown as horizontal dotted lines. The horizontal dashed 
line is the velocity at which the instability occurs in the GP 
simulation of the experimental 3D condensate of Ref. 0. 

order parameter develops a complex structure, starting 
from its center. This critical velocity is also shown in 
the lower part of Fig. [!|] (horizontal dashed line). It turns 
out to be very close to the one predicted for dynamical 
instability of the cylindrical condensate and well above 
the threshold of energetic instability. Finally, the critical 
velocity found in our GP simulations nicely agrees with 
the velocity at which the breakdown of superfluidity was 
observed in Ref. 0] (~ 0.5vb), therefore confirming that 
the dynamical instability plays a major role in this kind 
of experiments, as previously argued by Wu and Niu Q • 
The GP simulation also reproduces the shape of the 
density distribution after the onset of the instability (see 
insets of Fig. an d the density plots in Fig. 111(1 which 
is characterized by a condensed/coherent part plus a 
broader incoherent distribution on the right side, as ob- 
served in the experiment 0. The instability starts close 
to the center of the condensate. Radial excitations are 
involved in this process, as can be deduced by the occur- 
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FIG. 10: Center-of-mass velocity of the condensate during the 
its evolution in the harmonic trap after an initial displacement 
of 60 /im. The velocity initially increases, as expected for the 
harmonic dipole oscillation of the whole condensate, up to a 
critical value at which the condensate breaks up. The cor- 
responding quasimomentum p lies in the dynamical unstable 
region of Fig. In the insets we plot the linear density of the 
condensate along z at different times. The horizontal ruler 
has a width of 280 /im and the center of the trap is at mid 
point. 
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FIG. 11: Density plots in the (z,x) plane for different evo- 
lution times, corresponding to the inset pictures of Fig. 1101 
(t = 0, 30, 35, 40, 50 ms), after an initial displacement 
Az = 60 /jm. The parameters are those of the experiment in 
Ref. 0. 



rence of density patterns with radial nodes. 

By repeating the same simulations for much larger ini- 
tial displacements we find that the system can enter well 
inside the dynamically unstable region before decoher- 
ence processes take place. In this case, the acceleration 
of the condensate is so fast that dynamically unstable 
modes have not enough time to grow in a significant way. 
Eventually, for very large displacements the condensate 
passes through the entire dynamically unstable region 
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and reaches the zone boundary, at p ~ qs- At this point, 
the whole condensate undergoes a Bragg reflection soon 
followed by a sudden disruption process, involving soli- 
tons and vortex rings. This agrees with the recent GP 
simulations of Ref. 23], where the value Az = 150 
has been used and the nature of this instability has been 
discussed in detail. 



VI. CONCLUSIONS 

A general discussion of energetic and dynamical insta- 
bilities of a Bose-Einstein condensate moving in a ID op- 
tical lattice has been presented using the Gross-Pitaevskii 
theory. We have shown that transverse excitations are 
important in characterizing the stability diagram and 
the occurrence of a complex radial dynamics above the 
threshold for dynamical instability. 

The results of the full 3D calculations have been com- 
pared with those of an effective ID model, the nonpoly- 
nomial Schrodinger equation (NPSE), which includes the 
effects of the transverse direction through a Gaussian 
ansatz for the radial component of the order parameter 
[l5j . This model is shown to give accurate predictions for 
the instability thresholds, that are mainly determined by 
the dispersion of the lowest branch of excitations, with 
no radial nodes. 

The linear response analysis, in combination with the 



direct solution of the time-dependent GP equation, pro- 
vides a realistic framework to discuss the dissipative dy- 
namics observed in recent experiments. As an example 
here we have considered the case of experiments of Ref. 
0, providing a convincing evidence that the observed 
breakdown of superfluid flow is associated with the onset 
of a dynamical instability. 

The present analysis is also relevant in connection with 
the more recent experiments of Ref. |24| where the con- 
densate is loaded adiabatically in a moving lattice. In 
this case, interesting results have been found by com- 
paring the NPSE results with the experimental data for 
the nontrivial dissipative behavior of condensates loaded 
in higher Bloch bands |24| . In this respect, our work 
provides a further useful information, by pointing out in 
which cases, and for which observables, the NPSE is a 
reliable approximation of the full GP theory. 
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